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Abstract. Finite Element mesh generation remains an important is- 
sue for patient specific biomechanical modeling. While some techniques 
make automatic mesh generation possible, in most cases, manual mesh 
generation is preferred for better control over the sub-domain represen- 
tation, element type, layout and refinement that it provides. Yet, this 
option is time consuming and not suited for intraoperative situations 
where model generation and computation time is critical. To overcome 
this problem we propose a fast and automatic mesh generation technique 
based on the elastic registration of a generic mesh to the specific target 
organ in conjunction with element regularity and quality correction. This 
Mesh-Match-and-Repair (MMRep) approach combines control over the 
mesh structure along with fast and robust meshing capabilities, even 
in situations where only partial organ geometry is available. The tech- 
nique was successfully tested on a database of 5 pre-operatively acquired 
complete femora CT scans, 5 femoral heads partially digitized at intra- 
operative stage, and 50 CT volumes of patients' heads. In the latter case, 
both skin and bone surfaces were taken into account by the mesh reg- 
istration process in order to model the face muscles and fat layers. The 
MMRep algorithm succeeded in all 60 cases, yielding for each patient 
a hex-dominant. Atlas based. Finite Element mesh with submillimetric 
surface representation accuracy, directly exploitable within a commercial 
FE software. 

Keywords: mesh generation, Finite Element Method, elastic registration, 
mesh repair. 

1 Introduction 



Physically based models are now widely used in the field of biomedical en- 
gineering to represent human organs' geometrical and mechanical behaviors. 
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Among them, mimerical models based on the Finite Element Method [T] be- 
came very popular because of their ability to address the complex geometries, 
the anisotropic material properties and the specific boundary conditions asso- 
ciated with living tissues. 

In the field of medical imaging, the first Finite Element (FE) models were 
mostly used to better understand and validate a given surgical treatment, to 
model physiological behaviours or to provide virtual simulators for clinicians. In 
these frameworks, models were limited to a single generic model for each study 
(the terminology "Atlas" was often used). 

More recently, applications in the domain of Computer Assisted Planning 
and Computer Aided Surgery sparked the need for patient-specific FE models 
representing the modeled organ geometry reconstructed from patient medical 
image data, such as computed tomography (CT) or magnetic resonance imaging 
(MRI). In most cases, organs are identified in these data sets by means of manual, 
semi-automatic or automatic segmentation tools that extract shape information 
(3D points, contours and/or surfaces) necessary for the generation of the Finite 
Element mesh representing the volume of the organ. 

For both segmentation and FE mesh generation phases, manual intervention 
is often required which can make this procedure long and tedious. This is espe- 
cially true in situations where the pre- or intraoperative time window or clinician 
availability to perform these delicate tasks is limited. 

This paper addresses the second phase, namely the FE mesh generation, 
with the introduction of a procedure - the Mesh-Match-and-Repair (MMRep) 
algorithm - that allows a fast and fully-automatic patient-specific FE mesh gen- 
eration. 

Although FE models were already used in the field of biomedical engineering 
at the beginning of the 1980s [213] . researchers only began focusing on automatic 
mesh generation in the early 2005 [4i5.6.7„8ii9ll0lllll2j . The primary motivation 
for automatic mesh generation algorithms was that patient specific FE mod- 
els could be routinely used by the clinicians. In that domain, the vast major- 
ity of automatic mesh generators for living tissues produce tetrahedral meshes 
|13I14I15I16| . Based on a 3D surface representing the external geometry of the 
organ, these algorithms produce a volume made of high quality 3D tetrahedra 
and sometimes allow for adaptive mesh refinement ([E], Tetgero). 

Some researchers in biomechanics continue to propose hand-made FE meshes 
|18I19I20] . mainly for two reasons. Firstly, they argue that it is important to be 
able to identify sub-regions associated to anatomical sub-structures inside the 
3D FE mesh (e.g. the ventricles, the tumor and the hemispheres of the brain FE 
mesh proposed by Wittek et al. in ^20"). These sub- regions, corresponding to sets 
of elements, are labeled and associated with specific constitutive behaviors and 
boundary conditions. Secondly, they tend to prefer hexahedra over tetrahedra, 
based on numerical considerations [21] as well as the fact that for incompressible 
and/or nearly incompressible materials, 4-noded tetrahedra with linear shape 
functions tend to lock and become overly stiff [1]. Despite some improvements 

* |http: //tetgen.berlios .del 
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to this method of creating FE meshes this work stiU requires an excessive amount 
of manual effort to achieve satisfactory resuhs. 

In order to still benefit from this manual design while providing automatic 
FE mesh generation, techniques termed "registration methods" or "morphing 
methods" were recently proposed [4|22|11|12] . The main premise is to start with 
a predefined "generic" (or "template" ) FE mesh that represents the organ. This 
mesh is manually designed to include any necessary sub-regions and hexahedra, 
and to preserve element quality, orientation and density in regions that require it 
for numerical simulation. This template is then automatically morphed onto the 
patient data (3D landmark points, contours and/or surfaces) that was extracted 
from the segmentation of the medical images. This process generates a patient- 
specific mesh adapted to the geometries of the anatomical structures extracted 
from patient data, with a mesh topology that is similar to that of the template 
(same nodes and elements organization). 

Our group was the first to initiate this principle of mesh morphing with the 
introduction of the Mesh-Matching algorithm [4; . Since then, we encountered 
two strong limitations with our mesh morphing principle. Both are due to the 
fact that the template mesh quality after registration can be strongly decreased 
when the morphing algorithm induces excessive spatial distortions. Therefore, 
the first consequence is the inability to maintain the regularity of some elements, 
which disables FE analysis from being carried out. This concern was partially 
discussed in The second consequence of the morphing method is that the 
elements shape qualities are decreased in some regions of the template mesh, 
which leads to lower accuracy in the numerical simulation [23124125] . 

This paper aims at introducing our latest algorithms concerning (1) the elas- 
tic registration method that guarantees a C^-diffeomorphic transform; and (2) 
the mesh repair technique that ensures that the produced mesh complies with 
both regularity and quality criteria. The MMRep approach was applied and eval- 
uated on 60 clinical cases, which is, to our knowledge, the largest database ever 
tested in the literature for a patient-specific FE mesh registration method. 

The MMRep algorithm is a two-step sequential procedure. Firstly, the pa- 
tient data and the Atlas mesh surface nodes are registered using the elastic 
deformation procedure described in ^J2j This deformation is then applied to 
the inner Atlas nodes yielding a FE mesh that represents the modeled domain 
with sufficient accuracy. As a consequence of this deformation, the Atlas elements 
may suffer distortions and become either "irregular elements" which make FE 
analysis impossible, or "poor quality elements" in which case the computation, 
although feasible, can exhibit numerical instabilities. To recover mesh regularity 
and reach an acceptable quality level, the mesh repair procedure described 
in ^ is carried out on the deformed mesh. 
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2 Elastic registration 
2.1 Registration overview 

We define an elastic registration function as a mapping i? : M'' — > R'^ that super- 
imposes a source point cloud S onto a target, or "destination", data set D, which 
can either be a point cloud or a surface mesh. The computed elastic registra- 
tion procedure complies with continuum mechanics conditions on motion [26] as 
R defines a C^-diffeomorphic, non-folding and one-to-one correspondence 
between geometric data sets, as demonstrated respectively in ii2.3l and H2.4\ 

The input source points set is initially embedded in a deformable "virtual 
elastic grid" . We arbitrarily set the shape of the grid to be the bounding box 
of the points, extended by a 10% margin. The considered deformation R is 
formed by successive elementary grid deformations noted r, all having the desired 
regularity properties, much in the same way as proposed in [27'. 

The regular grid is progressively refined in order to increase registration ac- 
curacy and the expression of the compound registration function is: 

R = Vj' o . . . o r J o . . . o o . . . o 

where Gj,j = 1, . . . , J are successive grid refinements, and r*, i = 1,. . . ,Nj, 
elementary deformations performed at grid level j. As level indications are ir- 
relevant in the following demonstrations, the deformation R is simply written 
as: 

R — rjy o . . . o n (1) 

where N = J2j=i ^j- 

If required, the inverse registration, R~^, can be computed by combining the 
elementary inverses in the reverse order of the direct registration, thus: R~^ = 
rf ^ o . . . o TjY^ 

At each step n, the choice of the elementary deformation r„ to be applied 
to the source data is driven by the minimization of a "registration energy" E 
which measures the similarity between the deformed source points set and the 
destination data set D. As geometrical shape similarity is sought, E is defined 
as the sum of Euclidean distances between S and D: 

o i?„_i ) =^d(r„(i?„_i(s)),7?) (2) 
ses 

where = r^-i o . . . o ri represents the registration function assembled 

at step n — 1, and the operator •) can either be a point-to-point or a point- 
to-surface distance measure, depending on the nature of D. 

In order to speed up the computations of energy i?, a distance map is gen- 
erated from the destination data set D prior to registration [28]. Distance map 
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voxel dimensions are 1 x 1 x Imin in order to achieve submillimetric surface rep- 
resentation. The sampled space region covers the bounding box of the considered 
destination data, extended by a 5% margin. 

Point-to-surface registrations are performed in all three use cases discussed 
in §31 For each destination Atlas mesh, signed point-to-surface distances are 
measured using surface orientation information. A positive distance is recorded 
in the distance map for points lying outside the closed destination surface, and 
a negative distance is recorded for points lying inside. The distance map com- 
putations can take several minutes, which is not a limitation to our approach, 
even in the intraoperative context illustrated in §4.21 as the Atlas meshes are 
processed pre-operatively. At registration time, source-to-destination distance 
evaluations are done by trilinear interpolation within the signed distance map 
and the absolute value of the result is retained. 

The r„ functions are successively chosen so that E decreases optimally at 
each step. To this end the virtual grid is subdivided into a number of regular 
hexahedrons called "cells" . Each node of the grid is considered separately and the 
gradient of the registration energy is computed as function of the node's position. 
The opposite of this vector defines the node's preferred displacement, that is, 
the one that will lead to the greatest registration energy decrease achievable by 
moving the considered node. 

Let So :— S he the initial source points set, and let Si be the source points 
set at iteration i. Of all nodes, the one leading to the highest energy decrease is 
chosen and its preferred displacement is applied while all the other nodes remain 
fixed. This nodal displacement is propagated throughout the neighboring grid 
cells and the affected source points are moved accordingly to generate Si+i, the 
new set of deformed source points. 

Once this basis deformation step applied, the virtual grid returns to its initial 
regular configuration and the source points set Si+i is embedded within it. A 
new iteration can be computed by taking the new configuration Si+i as input. 
This Eulerian strategy allows large deformations of S to be achieved without 
having to maintain large-strain consistency of the virtual mesh, as would be the 
case if the virtual grid strictly followed the deformation with each iteration. 

The regularity of the grid before each iteration also saves computational time 
by allowing the interpolation of a node displacement throughout its neighboring 
cells to be computed using a template unitary cell. The iterations stop when 
no significant energy decrease can be achieved by moving any grid nodes at 
the current grid refinement level. In our implementation, we have set this stop 
threshold Te to 1%, and the registration iterations continue as long as energy E 
can be reduced by more than 1%, i.e. {Ei — Ei+i)/Ei > Te- Thus, the number 
of iterations performed at each refinement step is variable. 

The above iterative loop describes the procedure at a given grid refinement 
level. In order to maintain the spatial consistency of S during the assembly 
of the registration transformation, a top-down hierarchical approach has been 
implemented. The iterative assembly of the registration function R starts at the 
coarsest grid level. Once the deformation search at the current level has been 
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Fig. 1. Elastic registration overview, (a) at refinement level 1. (b) S'l after 
deformation at level 1. (c) Si at refinement level 2. (d) ^2 after deformation at 
level 2. 



exhausted, the grid is refined by subdividing each cell into 8 smaller ones in an 
octree method. 

Figure [T] illustrates in 2D the multi-grid iterative registration technique. The 
source points S are represented by the grey dots; for clarity D is not shown. In 
(a) the initial set So := Sis embedded within the virtual square grid discretized 
at the coarsest level 1, and the energy gradients are computed at the 4 mesh 
nodes. In (b) the optimal node displacement is applied to the grid and the source 
points move producing a new source set Si . If no significant energy decrease can 
be generated at level 1, the grid is refined at level 2, (c), 5"! is embedded within 
it, and the energy gradients are computed at the 9 nodes of the mesh. In (d) 
the best nodal displacement is applied and the resulting source points set S2 is 
computed. 

At each grid refinement level, the size of the registered source features is 
approximately that of the current grid cell. The virtual grid nodal displacement 
is only applied if the resulting grid deformation leads to a registration energy 
decrease. As all the source points located in the grid cells surrounding the dis- 
placed node are altered, the source features significantly smaller than the current 
cell size have limited impact on the registration sequence. The dimensions of the 
smallest cell reached during the top-down refinement descent thus roughly define 
the size of the specific features present in S that may not be registered if the 
corresponding feature is absent from D. 

By primarily focusing on larger source features, this hierarchical approach 
also reduces the influence of noise usually found in patient data gathered from 
pre-operative medical image segmentation (see i 34.1|) or from intraoperative ac- 
quisitions (see M.2|) . During the evaluation of the MMRep technique the smallest 
cell size was set to 1mm, which required about 8 or 9 grid subdivision levels as 
the typical model size was about 25cm and 250mm/2^ « 1mm. 

Furthermore, the mechanical regularization strategy described in >j2. 61 limits 
excessive space distortions due to the presence of noise by monitoring the magni- 
tude of potential elastic energy in the source space throughout the deformation 
process. Accurate registration of outliers has indeed a prohibitive mechanical 
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cost which, in our procedure, discards elementary deformations attempting to 
register them. 

The optimization procedure is done by a gradient descent technique 1291 . The 
energy gradients are evaluated at each node by the Finite Differences method. 
The line search in the direction of the gradient descent is performed by approx- 
imating the energy curve by a parabola P(t) — at^ + ht -\- c. The parameters 
a, b and c are deduced from: the current value of E at the considered node, 
corresponding to P(0); the slope of E in the descent direction, defining P'(0); 
and the value of E at the maximally displaced node position (see yield- 
ing P{1). Finally, the optimal descent step in the current direction is given by 
tmin = -b/{2a). 



2.2 C^-difFerentiability of the deformation 

At each registration step, the displacement applied to a given grid node n is 
propagated to the source points located in cells surrounding n by means of a 
"weight" function w^:B? ^ [0, 1]. 

Let it be the displacement applied to node n, and s be a source point found 
in a cell affected by the movement of n. The displacement propagated to s is 
Wn{s,)lt . Similarly to the shape functions in the Finite Element Method, the sup- 
port of the weight function is the union of the cells neighboring n. The defor- 
mation consistency is further ensured by the two following conditions: u'n(n) = 1 
and w„ = at the boundary of its support. 

From the above, it follows that the elementary elastic registration function r 
created by the displacement it of node n has the following expression: 



3 ^ R3 



, S l-^- S + Wn(s)lt (3) 

C^-differentiability of the elementary registration function r stems from 
the C^-differentiability of the weight function Wn and the Jacobian matrix of r 
is given by: 

Now, let So G M"^ be an arbitrary source point and s„ this point transformed 
after applying n elementary registration steps. Then {s„, n G [0, N]} is the set of 
all successive positions of Sq and, according to Eq. [U R{sq) = s^r. The Jacobian 
matrix of R computed at point Sq is given by the following chain rule: 

Jr{so) := ^(So) = ^— (SAr_i) ... -^(si) ^(sq) (5) 

as as as as 

All weight functions stem from a "template" weight function "w" defined 
on the [0, 1]'^ grid cell and associated to node (1, 1, 1). Let tt be a third degree 
polynom defined by 7r(f) = t'^{3 - 2t), such as 7r(0) = 0, 7r(l) = 1, 7r'(0) = and 
7r'(l) = 0. The template weight function w is defined on [0, 1]'^ as: 
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Fig. 2. (a) Variable changes required in 2D for the assembly of Wn on a neigh- 
borhood of 4 cells centered on node n = (1, 1). (b) Value of Wn plotted over the 
four 2D cells, (c) An elementary deformation leaves all segments [a, b] parallel 
to the applied nodal displacement unchanged. 



w{s) = w{si, 82,83) := 7r(si)7r(s2)7r(s3) (6) 

A specific weight function is defined on the union of neighboring cells 
around the displaced node n by variable change and scaling in order to adapt the 
canonic [0, 1]'^ domain to the cell size within the actual grid. Fig. [5]-a illustrates 
on a 2D grid the variable changes needed to construct from w the weight function 
Wn, with n = (1, 1), defined on a 2x2 cells neighborhood of n. Fig. [2]-b shows 
the 3D plot of the resulting weight function Wn ■ [O7 !]■ 

The two other regularity properties of the elementary registration functions 
r, namely bijection and non-folding, are enforced by limiting the amplitude 
of nodal displacements, as described in the following sections. 

2.3 Control over space distortion 

Space folding can be mathematically expressed as follows. Consider a positively 
oriented set of three infinitesimal vectors, dXi, dX.2 and 1^X3, placed in unde- 
formed space at point s and defining the volume dV. The signed value of dV 
can be computed as the determinant of the matrix formed by the three vectors. 

dV = {dXi X dXa) • dXs = | dXi dXa dXa | 

R is said to be locally non-folding if the deformed infinitesimal volume 
dv = R{dV) defined by the three deformed vectors dxi, dx2 and dx3, remains 
positive, or, in other words, if R does not change the orientation of space in the 
neighborhood of s. 

The differential form of R gives us the expression of the deformed vectors: 



Vi = 1,2,3, dx, = — (s) dX, 
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and the relation between dv and dV is: 



dv 



I dxi (ix2 (ix3 



OR 



IdXi dX2 dXgl = \Jr{s)\ dV 



where J_r(^) •— {dR/ds){s) is the Jacobian of R given by the chain rule in 
Eq. [SJ It follows that an application R is non-folding if its Jacobian is strictly 
positive, i.e.: 



Vs e r3, \ Jti{s)\ > 



(7) 



If the value of the Jacobian at a given point is greater than 1, the deforma- 
tion locally stretches space and if its value is smaller than 1, space is locally 
compressed. 

From Eq. [5]it follows that R is non- folding if and only if all r, , i = 1, . . . , 
are non-folding, since: 



\Jr\ 



dR 
ds 



N 

n 



drj 
ds 



N 



Let's now see how to ensure that each individual elementary registration 
function is non-folding. Using Eq. SI the non-folding condition on r becomes: 



Vs e |./^(s)| = 1 -I- ^w;„(s) • ^ > 



(8) 



Given the above expression of the template weight function w, the magnitude 
of the gradient of a specific weight function w„ defined on a [0,L]'^ grid cell has 
an upper bound of 3a/3/(2£). To ensure non- folding it is thus sufficient to limit 
the amplitude of nodal displacements: 

<38%L< 



3^3 



In our implementation of the algorithm this value has been reduced to 10% 
of the current grid cell size, so as to not only achieve non-folding but also 
limit space distortion at each elementary registration step. As < L/10 and 
ll^Wnll < 3V3/(2L), it follows that: 



3%/3 
"20" 



< 



l+^Wn(s)-lt < 1 + ^ 



(9) 



The above equation gives us, for all cell sizes, approximatclfl Jacobian lower 
and upper bounds: 0.74 < \Jr\ < 1.26, and hence 0.74^ < \Jr\ < 1.26^. 



^ 3^3/20 ~ 0.259 
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2.4 Registration inversion 

In this section we will show that under the non- folding constraint, the deforma- 
tion i? is a bijection and discuss how its inverse, R~^, can be computed with a 
pre-defined level of accuracy for any point q e M'^. 

The registration function R is one-to-one if the same is true of each elemen- 
tary registration r. To prove that a non-folding elementary registration function 
is also a bijection, consider a 2D deformation r, as defined in Eq. [31 created by 
applying the displacement it to the central node of the 4 cells group depicted 
in Fig. [2}c. 

Given that all the displacements applied to the source points are collinear 
with it, parallel segments such as [a, b] in Fig. [21-c, map onto themselves. Fur- 
thermore, points lying outside or on the boundary of the 4-cell group are left 
unchanged. To prove that r is a bijection it is thus sufficient to show that it 
defines a one-to-one application [a, b] — ^ [a, b] , for any segment [a, b] parallel to 
It. 

As the considered segment is parallel to it, there exists a scalar /3 such 
as b = a -I- f3lt . Moreover, all segment points p S [a, b] can be written as 
p = a -I- pit, with p G [0,(3], leading to r(p) being rewritten as follows: 

''(p) = P + Wn(p)^ = a + (p + Wn(a + plt))lt 

The deformation r is thus a bijection [a, b] — > [a, b] if and only if the mapping 
/ : p I-)- p -I- Wn(a + pit) is also a bijection [0, /3] [0, which we shall prove 
is true under the non-folding hypothesis. 

Deformation r leaves the segment extremities unchanged, r(a) = a and 
r(b) = b, which, using the above expressions, leads to /(O) = and = /3, 
respectively. Furthermore, the derivative of / is /' = 1 -|- ~^Wn - it . As a con- 
sequence, if r is non-folding then, according to Eq. [H we have /' > which 
implies that the strictly monotonic function / is a bijection [0, /3] — > [0,/3] and 
so is r : [a, b] — )• [a, b] , which concludes our proof. 

Another issue is computing the inverse registration — r^^ o . . .o with 
a user-defined accuracy e, measured in undeformed space. The task consists in 
finding, for a given deformed point q e M"^, an undeformed point p £ M"^, such 
as ||p — i?^^(q)|| < e. We will first show how to accurately compute the inverse 
of an elementary registration function r, and then use these results to accurately 
inverse the compound registration R. 

Given q G [a, b], q = st.-\- pit , finding p = a+ pit such as p = ?'^^(q) reduces 
to solving p = f~^{p) on [0,/3]. The mapping / is a 9**^ degree polynomial (see 
Eq. [6]) and as no analytical expression of the solution is available, it must be 
approximated iteratively using, for example, a Newton-Raphson procedure. 

Let p = a -f 'pit = r~^(q) be the searched point, and p = a -)- pit r~^(q) 
its iteratively computed approximation. The approximation error in undeformed 
space can be rewritten as: 



l|p-r-i(q)|| = ||p-p|l = |p-p|||lt|| 



(10) 
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In order to compute the inverse registration with the desired accuracy level 
e, the exact solution p = f^^(p) must be approximated by p so that: 



As the value of p is unknown, the approximation error must be computed 
in deformed space. To this end the above expression is rewritten and, using the 
finite increments theorem, yields: 

\p-p\ = \r\fip)) r\p)\ < M \f{p) p\ (11) 

The above constant AI is computed using the relation (/^^)' = {f')^^, which 
in conjunction with Eq.H gives 1/1.26 < (/"^)' < 1/0.74 = M . This leads us to 
the conclusion that in order to compute an e -accurate inverse of q in undeformed 
space, Newton-Raphson iterations must be carried out until the approximate 
parameter p satisfies: 

<0.74^ 

Now let's compute the inverse of a compound registration function R. To 
control the accumulated error overhead at each elementary inversion step we 
will use the fact, demonstrated by combining Eq. [TU] and [11] that the inverse of 
any elementary registration function r is Ai'-Lipschitz continuous, i.e.: 

Vq,p e R3, \\r-\q) - r-\p)\\ < M ||q - p|| 

For the sake of clarity, let R — o r2 o ri be the considered registration 
function and S3 G K'^ a point in deformed space which inverse, R~^{s^), is 
sought with accuracy e. Adopting the notation convention used to derive Eq. [5] 
above, we have: 



i?-i(s3) = (rfi o r^^ o r3-i)(s3) = (rf^ o r^^){s2) = r^\si) = Sq 

Fig. [3] illustrates the 3-step computation of i?^^(s3) along with error accu- 
mulation due to approximations performed at each elementary inversion step. 
Undeformed space is represented on the left, and deformed space on the right of 
the figure. Approximation steps are shown as dashed oblique arrows, and exact 
inversions as horizontal arrows. 

We shall now establish the relation between e, the selected tolerance in un- 
deformed space, and successive elementary approximation errors €3, £2 and ei. 

Inversion of r^. S2 = r3"^(s3) is approximated by t2 with accuracy es in 
r3-undeformed space: £3 > ||t2 — S2||. Using the Lipschitz constant M, this error 
propagates to the left as: £3 M > ||ti — Si||, and £3 A'P > ||to — So||. 

Inversion of 7'2- ti = r^^(t2) is approximated by Ui with accuracy £2 in r2- 
undeformed space: £2 > ||ui — ti || . As in the previous step, this error propagates 
to the left as: £2 M > ||uo - to||. 
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4^ 4^ ^ 

Fig. 3. Approximation of Sq = i?^^(s3) by Vq. Horizontal arrows represent ex- 
act inversions and oblique dashed arrows iterative inverse approximations. Left: 
undeformed space, right: deformed space. 



Inversion of ri. Finally, Uq = (ui) is approximated by vq with accuracy 
€i in initial undeformed space: ei > ||vo — Uo||. 

The final error in undeformed space ||vo — So|| can be upper bound as: 

||vo - Soil < ||vo - uoll + ||uo - toll + ||to - Soil <ei+€2 M + eg 

The above expression can be generalized to any compound registration func- 
tion R = Tat o . . . o n. In order to meet the desired accuracy standard e, the 
computational effort can be spread among all N elementary inversions by set- 
ting the Newton- Raphson approximation threshold for each ,n — I, . . . , N 
to: 

e 

This is a reasonable heuristic as e„ is smaller for higher values of n, cor- 
responding to finer grid levels and hence, smaller nodal displacements (10% of 
grid cell size - see t ^2.3p . Solution search intervals are thus narrower, and higher 
inversion accuracy is thus easier to achieve than at coarser grid levels. 

There is a final consideration on registration inversion: if the cubic polynomi- 
als TT used to define the shape function w in Eq. |6] are replaced with the identity 
i : [0, 1] — >■ [0, l],t t, then / becomes a 3'''* degree polynomial and the exact 
inverse of the registration function R can be computed in a single iteration by 
using the well-known Tartaglia-Cardan formula. Nevertheless, this enhancement 
comes at a price as the smoothness of the registration function R is degraded 
from Ci to C°. 



13 



2.5 Asymmetric registration 

The "asymmetric" energy E defined in Eq. [2] liandles situations wtiere ttie fea- 
tures present in S are present in D but tlie reciprocal is not necessarily true. 
A situation can occur in which the Atlas mesh only represents a fraction of the 
organ features available within the patient data (see M.'6\i . requiring that the 
elastic registration R is computed as transforming the Atlas to fit the patient 
data. The patient specific mesh is then obtained by application of this transform 
to the generic mesh i.e. Patient mesh = i? (Atlas). 

Conversely, if only partial information about the patient's organ is available, 
its overall shape has to be inferred from the a priori knowledge carried by the 
generic Atlas mesh (see H4.ll and 14.21) . This is done by taking the patient data 
as source points, and computing the registration function R that fits the source 
points onto the Atlas. The patient specific mesh is obtained by inverting the 
resulting elastic deformation, i.e. Patient mesh = i?"-*- (Atlas). 

2.6 Mechanical regularization 

In order to avoid excessive space distortions, a mechanical regularization term is 
monitored during the registration process, thereby upgrading the virtual elastic 
"grid" concept to virtual elastic "solid" . As the elastic registration compensates 
for inter-individual morphological variations and does not model a physical de- 
formation, the underlying ad-hoc mechanical properties are not related to the 
actual rheology of the organ under study. 

Using the notation from Eq. [21 we now define the Jacobian matrix of the 
overall registration function considered at iteration n, and taken at material 
point X, as: 

Equation [5] shows that this matrix can be updated after each elementary 
registration step by first multiplying the previous matrix with the Jacobian 
matrix of the new elementary deformation r„_|_i taken at the current position of 
the considered material point, (r„ o . . . o ri){X). 

During the registration assembly, the elastic energy stored in the virtual 
solid is measured at a set of material points evenly distributed among 

the initial source data, and set to "probe" the space distortions induced by the 
accumulated elementary deformations. To do so, at each iteration n the Green- 
Lagrange strain tensor e„ is derived from the above mentioned Jacobian matrix, 
as e„ = (Jj J„ — I)/2, and the stress tensor (T„ is related to the strain tensor e„ 
by a linear constitutive equation cr„ ij = Dij^i e„ ■ 

The potential elastic energy generated by the deformation iZ„ at control point 
A* can be computed using the Total Lagrangian formulation, as: 
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where is the volume in the initial source configuration associated to, or 
"monitored" by the material point X^. In order to preserve fast registration 
computation, the above integral is approximated as w |y*| e„(X') : cr„(X*), 
where \ V^\ is the measure in the initial configuration of the volume associated 
to the material control point X'. 

The sum of the contributions of all control points {X^}i gives an approxi- 
mation of the total potential elastic energy stored in the deformed source space 
at iteration n, as W„ = J2i^n- This measure is taken into account at each 
deformation step to select, among all possible elementary deformations, the re- 
gistration function r„ which offers the best ratio between registration energy 
decrease and elastic energy increase. 

To this end, at each iteration n and for each candidate deformations r„ 
the associated registration energy decrease AEn > is computed as AEn = 
E{Rn-i) — E{rn o Rn-i)- The change in potential elastic energy AWn associated 
to each r„ is also computed as AWn = Wn — W„_i and the following selection 
algorithm is applied. 

1. Among all candidate elementary deformations, only the deformations r„ 
leading to a relative registration gain greater than the stop threshold Te are 
considered, i.e. the ones for which AEn/En-i > Te- 

2. Among those, if deformations such as AWn < can be found, the one 
yielding the highest registration energy decrease is chosen. 

3. Otherwise, the deformation r„ having the highest AEn/ AWn ratio is applied. 

Rule 1 merely implements the stop criterion mentioned above. If no elemen- 
tary registration function satisfying this condition can be found, the iterations 
stop. Rule 2 favors space "decompression" if it goes with satisfactory registration 
enhancement. Finally Rule 3 is applied in the most general case to select the 
deformation offering the best trade-off between space distortion and registration 
gain. 

It is important to emphasize here that although the elastic grid initial state 
is restored prior to each elementary registration step (Eulerian approach), the 
mechanical energy terms W„ computed above keep track of all the deforma- 
tions accumulated in i?„ at step n (Lagrangian formulation). They are therefore 
appropriate for measuring space distortion as the registration advances. 

As a relationship between inter-individual shape variations and mechanical 
behavior could not be determined, the simple St. Vcnant-Kirchoff mechanical 
framework was selected. An isotropic soft compressible material using an empir- 
ically set Poisson's ratio u = 0.2 and a Young's modulus E = IPa was imple- 
mented. The value of E has no effect on the registration sequence as it is not the 
sum but the ratio between registration and elastic energy that conditions each 
elementary deformation. 

Finally, in our implementation, 1000 control points are used. The {X*}i=i....^iooo 
are distributed in the bounding box of the source data, on a regular 10 x 10 x 10 
grid and all \ V^\ are set to 1/1000 of the bounding box volume. 
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2.7 Multiple structures registration 

Prior to mesh registration, a subset of the Atlas mesh nodes needs to be labeled 
as anatomical features that can be identified within the patient data. 

In the case of the femur models (see i|4.1l and l4.2[) the surface mesh nodes are 
labeled as "bone surface" as they lie on the cortical surface of the bone, which 
can easily be extracted from a Computed Tomography (CT) volume. 

As for the face model, two families of nodes are defined, labeled "skin" for 
the exterior surface nodes registered onto the patients' faces skin, and "bone" for 
the interior surface nodes that correspond to the skull features and need to be 
registered onto segmented patients' skulls (see ^4^. This application example 
illustrates the capability of our procedure to capture the shapes of multiple 
anatomical structures modeled within a unique generic FE mesh. 

When multiple anatomical structures need to be recovered, the registration 
of the FE mesh onto the patient data is driven by the minimization of an energy 
term E^^^ which measures the fit between the labeled nodes and their corre- 
sponding anatomical structures: 

L 

1=1 seSi 

where I = 1, . . . ,L are the predefined labels and Si and Di are the corre- 
sponding source and destination regions respectively. When multiple labels are 
defined, a distance map is computed for each Di subset of the destination data. 

The elastic grid deformation is driven by the E^'^^ energy computed over 
the sets of labeled source nodes and destination surfaces. The unlabeled source 
points, on the other hand, are embedded within the elastic grid and passively 
follow the deformation without contributing to the energy term. 

3 Mesh repair 

Although the elastic registration algorithm described above strongly limits space 
distortions, the registered mesh may exhibit irregular or low quality elements 
that need to be untangled before proceeding to FE analysis. Indeed, the non- 
folding nature of the registration function is a local property ensuring that space 
orientation is preserved. While this is locally true at every point in space, the 
property does not hold when considering finite structures such as element edges 
and, as a consequence, irregular elements may appear after the application of the 
smooth and non-folding elastic deformation to the initially regular Atlas mesh. 

Element regularity is assessed by considering the mapping F : ^ i— ^ x, 
between the element parent (or reference) coordinates system (^1,^2,^3) and the 
actual element coordinates (a:i, X2, xa). The Jacobian J(^) is the determinant of 
the matrix dF/d^. Its value represents, at a given reference point ^, the local 
volume transformation between the parent and actual element configuration. 

An element is said to be regular if J(^) > for all ^ in the parent config- 
uration. Element regularity is usually assessed by considering the value of J at 
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specific points such as the integration points or the element nodes 1 . We call a 
node n "irregular" if it has a negative Jacobian J^, computed in element e. A 
mesh is said "regular" if all the elements, and hence all nodes, are regular. An 
irregular mesh is not suitable for Finite Element analysis as the singularity of 
the mapping F leads to modeling inconsistency. 

Element quality, on the other hand, is a measure of the conformity of its 
shape, which reflects the evenness of the discretization of the modeled domain. 
There is a great variety of quality measures and their relevance is dependent on 
the considered element type and computations to be carried out [23125) . 

Given the fact that fast and robust tetrahedral discretization solutions are 
already available, such as the widely used TetGen software, or have been pro- 
posed in the literature [13130115] . we focus our work on hexhedral-dominant 
meshes and demonstrate our mesh generation technique using a popular quality 
measure, well suited for hexahedrons and wedges: the Jacobian ratio (JR) |3T| . 

The JR is defined for a given node n considered within element e. Its value is 
the ratio between the Jacobian at node n in e, and the maximal nodal Jacobian 
value in element e, J^ax — Maxmge{>/m}; thus: 



By comparing, the Jacobian value of each node with the maximal value within 
the considered element, the JR gives an indication of the contribution of a node 
to the overall element distortion, as opposed to local nodal distortion between 
the parent and actual configuration measured by the Jacobian alone. Element 
quality with respect to parent configuration is proportional to the JR value, 
ranging from to 1. 

Mesh repair is a two-fold process: first all the elements in the mesh are 
inspected and, if necessary, their nodes' positions are adapted so as to recover 
regularity; then a second relaxation procedure is carried out on the mesh if 
the achieved quality levels are unacceptable. Both steps can affect the nodes 
positions and alter the organ surface representation achieved after the elastic 
registration step. 

In the absence of a formal proof that an acceptable mesh configuration exists 
and can be found by the relaxation procedure in all situations, the aim of this 
study is to evaluate the performance of MMRep on a database of diverse organ 
shapes and clinical use cases. The results described in ^ 14.11 14.21 and 14.31 suggest 
that the smoothness of the elastic registration strongly limits spatial distortion 
making it possible to recover both mesh regularity and quality without inter- 
fering with the prior surface registration, by applying small displacements to a 
limited subset of the surface nodes. 

A large number of nodes in a mesh can make the repair procedure com- 
putationally prohibitive. This complexity can be greatly reduced using a local 
relaxation strategy i.e. grouping all irregular or poor quality nodes into "re- 
gions" defined in such a way that nodal corrections applied inside a region leave 
the outside mesh configuration unchanged. This local repair strategy makes it 
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possible to perform all relaxation procedures independently on each repair re- 
gion identified within the mesh. It also decreases significantly the computational 
complexity as the number of degrees of freedom to be considered in a region is 
usually small. 

The untangling of an irregular region R consists in finding a configuration 
where all the Jacobians in R, { Jjjjgfl, are positive. This nodal relaxation can be 
formulated as a maximization procedure driven by a "regularity energy" Ej^ . En 
is expressed as the sum of all Jacobians within R affected by a penalty function 
tpk, of strength controlled by an index fc, giving: 

where </?fc(<) = 1 — exp{—kt). As the parameter k is increased during the 
optimization process, the infiuence of negative values overbalances the positive 
ones thus favoring a solution where all Jacobians in the sum are positive. 

The aim of the quality maximization, on the other hand, is to find a con- 
figuration where all the JRs in a region R are above a predefined level JRmi„. 
The associated "quality energy" Eq is defined as the sum of the JRs within R 
weighted by a penalty function thus: 

jeR, e£R 

where ipkit) = 1 — exp{k{JKmin — t))- Similarly to the regularization process, 
the penalty parameter k is used to find a solution where all JRs contributing 
to the sum Eq are above the quality threshold JRmin • The value of JRmm was 
chosen in accordance with the quality standard requested by the commercial FE 
analysis software ANSYS Workbench (ANSYS Inc., USA) i.e. JR,„,„ = 1/30 
[32]. 

The initial value of k for both penalty functions ipk and ijjk must be care- 
fully chosen. Indeed, given the formulation of the penalty functions, an excessive 
penalization level may induce strong numerical instabilities in the optimization 
process. To avoid this issue, the starting value of k is determined by considering 
the slope of the penalty function at the most penalized energy term (minimal 
Jacobian during regularization phase; minimal JR during quality optimization 
phase) and ensuring that it does not exceed a predefined threshold. 

Both regularity and quality optimizations are carried out by gradient ascent. 
Gradients of both Er and Eq energies are computed using the centered differ- 
ences scheme. Assuming unimodality of the local energy function, the maximum 
search in the direction of ascent is done using the golden section technique [29] 
between the current nodal configuration and the configuration obtained after 
applying the maximal amplitude correction. 

The amplitude and number of iterations for each mesh region are limited so as 
to restrict loss of surface representation accuracy. Our implementation allows a 
maximum of 50 iterations, and nodal displacements at each step have a maximal 
amplitude of 0.1mm. The 50 x 0.1 = 5mm maximal displacement can only be 
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achieved by moving a unique node in a constant direction throughout the repair 
procedure, which seldom occurs. During the experimental validation discussed 
in 21 nodal corrections with mean amplitude less than 1.2mm and applied to 
less than 1% of the nodes were sufficient to repair the 60 registered meshes. 

The value of maximal nodal correction used here is not universal and must 
be determined for each field of application according to the dimensions of the 
modeled domain and maximal tolerance on the representation of its geometry. 
In our case, the results were consistent with the desired submillimetric mean 
surface representation accuracy, as shown in tables (TJ [3] and [5] 

4 Experimental evaluation 

The MMRep technique was evaluated in three different situations where: 

— complete organ geometry can be retrieved from the available data, as 
discussed in !j4.1l 

— only partial organ geometry is accessible, as illustrated in H4:.2\ 

— distinct organ features need to be taken into account by the generated 
model, as shown in 

The mesh adaptation technique alone is presented here and the discussion 
about the biomechanical simulations illustrating the three following "use cases" 
falls out of the scope of this article. 

4.1 Complete pre-operative femora CT scans 

In this section, we demonstrate our technique in the context of total knee arthro- 
plasty (TKA) . The prosthesis placement can be optimized to avoid unsealing or 
femur fracture using biomechanical modeling and FE analysis of the stresses 
within the tissues. These tissues are modeled based on the bone mechanical 
properties, geometry and prescribed loads inferred from the patient morphology 
and gait [55] . 

Mesh registration procedure The manually assembled right femur Atlas 
mesh [34] used in this part of the study is composed of 4052 nodes, forming 3018 
elements : 2960 hexahedrons and 58 wedges (6 nodes prisms). The elements are 
organized so as to reflect the bony structure such as the femoral diaphysis cortex 
which is discretized by a single layer of elements. The principal mesh features 
are illustrated in Fig. |4l 

For the 5 considered patients a CT scan of the right leg was pre-operatively 
acquired and a semi-automatic threshold-based segmentation procedure was per- 
formed by a clinician in order to identify the femur cortical surface. Then the 
resulting points cloud was rigidly registered onto the Atlas bone surface using 
an Iterative Closest Point algorithm [55] . 

Each elastic mesh registration was carried out as described in [J2] by taking 
the source points set 5 as Si, the segmented cortical points for patient i, and 
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(d) 



Fig. 4. Femur Atlas mesh from Couteau et al. (a) Overview, (b) Femur 
head and great trochanter, (c) Cut out showing the diaphysis cortex layer of 
elements, (d) Distal condyles. 



the destination set D as the Atlas mesh cortex. Once the segmented points were 
registered onto the surface of the Atlas, the application of the inverse registration 
function R^^ to the generic mesh produced the patient specific FE model. The 
chosen accuracy level for the inverse computation, as described in ^ ^2.41 was 
e=0.1mm. Mesh regularity and quality criteria were subsequently analyzed and 
the mesh repair procedures described in f}3]were applied to the model. 



Results Table [T] describes the performance of the mesh registration procedure 
for the 5 data sets. The surface representation error is computed as the distance 
between the segmented CT points and the generated patient specific FE mesh 
surface. The registration times include the direct registration R computation as 
well as the application of the inverse registration function R^^ to the Atlas mesh 
nodes. 

Table [5] gives the performance of the mesh repair procedure carried out on 
the 5 deformed meshes. Repair computation times and numbers of corrected 
nodes are given along with nodal displacements statistics. 

The MMRep procedure succeeded in generating a quality femur mesh for all 
5 patients. All computations were carried out in less than one minute, which is 
acceptable even in an intraoperative context. The surface representation figures 
remained unchanged after the application of the mesh repair procedure (see 
Table[T]) as the proportion of displaced nodes did not exceed 1% and most of them 
were inner nodes which could be freely moved without affecting the mesh surface 
shape. The reported maximal errors are mainly due to manual segmentation 
irregularities and lack of local refinement of the Atlas mesh, making it difficult 
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Patient 


1 


2 


3 


4 


5 


Registration (sec) 


25 


36 


42 


26 


32 


CT points 


10930 


25980 


22924 


20065 


17886 


Mean err. (mm) 


0.3 


0.4 


0.3 


0.4 


0.3 


Max err. (mm) 


5.4 


6.6 


5.2 


6.0 


5.5 


er (mm) 


0.5 


0.5 


0.4 


0.5 


0.4 



Table 1. Registration (sec): elastic registration times, in seconds; CT points: 
quantity of segmented points in CT volumes; Mean, Max, cr (mm): surface rep- 
resentation mean and maximal error, standard deviation, in millimeters. 



Patient 


1 


2 


3 


4 


5 


Regularity (sec) 


2.2 


0.6 


0.4 


0.9 


0.5 


Quality (sec) 


3.6 


2.6 


6.8 


0.9 


0.7 


% nodes 


1.0 


0.3 


0.6 


0.5 


0.4 


nodes/4052 


39 


27 


23 


22 


17 


Mean disp. (mm) 


1.0 


0.3 


0.2 


0.5 


0.2 


Max disp. (mm) 


3.4 


2.5 


1.1 


1.9 


0.6 


(T (mm) 


1.0 


0.5 


0.2 


0.5 


0.2 



Table 2. Regularity, Quality (sec): mesh repair times for both phases, in seconds; 
% nodes, nodes/4052: fraction and number of nodes moved by the repair pro- 
cedure; Mean, Max, a (mm): mean and maximal nodal displacements, standard 
deviation, in millimeters. 



to capture some local shape variations. If necessary, this issue could be solved 
by further refining the generic mesh. 

A sample result of the procedure is shown in Fig. [5] with focus on proximal 
and distal parts of the automatically generated femur mesh. 

Mesh quality distribution measured on the 3018-element full femur Atlas 
mesh is shown in Fig. |6]-a, and Fig. |6]-b gives the mean quality distribution in 
the 5 generated patient-specific meshes, along with standard deviations. 

The histograms presented in this article, in Figs. El |8] and [l2l classify the 
elements into 5 Jacobian ratio categories. The first interval [0.03,0.2] lists the 
elements with a quality level that is acceptable from the point of view of our 
target application ANSYS, but is usually considered to be "questionable" . Due to 
the manual assembly process, a small number of questionable elements is present 
in the Atlas meshes used in this study. As shown by the figures, this proportion 
only slightly increases after the application of the elastic deformation and mesh 
repair procedures. 

4.2 Partial intraoperative femora digitizations 

In this part of the study, we demonstrate the possibility of intraoperative FE 
mesh generation during a total hip arthroplasty (THA) procedure. As in TKA, 
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FE analysis can be used here to optimize femoral stem placement so as to min- 
imize internal stresses and maximize the implant lifetime p36] . 

Mesh generation procedure Due to the surgical procedure limitations, the 
complete shapes of the patients' femoral heads were not available pre-operatively 
and each bone geometry was acquired intraoperatively by sliding a calibrated 
pointer on the cortical surface of the partially exposed hip. The pointer position 
was tracked in space by means of an optical localization system (Polaris, NDI, 
Canada) and the position of its tip was recorded continuously. 

The initial positions of the recorded points with respect to the Atlas model 
were computed using the correspondence between anatomical landmarks such 
as the knee center and the piriformis fossa, localized in patient space using the 
pointer, and defined in the Atlas space by an expert operator. Fig. [7] shows the 
similar and very localized distributions of the digitized points for all 5 patients 
with respect to an approximative surface model of the proximal femur (3 right 
and 2 left hips were included in this study). 

The hip Atlas mesh used here is the upper part of the complete right femur 
Atlas presented in §4.11 After truncation, the right hip model was mirrored with 
respect to the sagittal plane to produce the left hip generic model. It is composed 
of 2105 nodes, forming 1738 hexahedrons and 16 wedges organized so as to reflect 
the hip principal mechanical structures such as the femoral head and neck. 

Mesh registration and repair was carried out as described in HA.ll taking 
S as Si, the digitized points cloud for patient i, and D as the Atlas hip bony 
surface. Each patient specific mesh was created by applying the inverse R-^ of 
the registration function computed with an accuracy level e=0.1mm, as described 
in ijMl 

Results The two result tables below are similar to those included in Ta- 
ble [3] gives the performance of the mesh registration procedures carried out on 
the 5 data sets, the surface representation error being, this time, computed by 
considering the distance between all digitized points and the generated patient 
specific hip surface. 



Patient 


1 


2 


3 


4 


5 


Registration (sec) 


15 


17 


20 


11 


23 


Points 


1204 


1433 


1421 


1450 


1437 


Mean err. (mm) 


0.3 


0.4 


0.4 


0.3 


0.3 


Max err. (mm) 


2.1 


3.4 


4.2 


2.2 


2.7 


(T (mm) 


0.3 


0.4 


0.4 


0.3 


0.3 



Table 3. Registration (sec): elastic registration times, in seconds; Points: number 
of intraoperatively digitized points; Mean, Max, a (mm): surface representation 
mean and maximal error, standard deviation, in millimeters. 
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Table |4] presents the performance of the mesh repair procedures for the 5 
meshes. The repair times are given in seconds. In the case of patient 4, the 
patient mesh produced by the elastic registration was already regular. 



Patient 


1 


2 


3 


4 


5 


Regularity (sec) 


0.4 


0.1 


0.5 





0.4 


Quality (sec) 


0.5 


0.4 


0.8 


0.4 


0.4 


% nodes 


0.3 


0.05 


0.7 


0.05 


0.4 


nodes/2105 


6 


1 


14 


1 


8 


Mean disp. (mm) 


1.2 


0.06 


0.6 


0.04 


0.7 


Max disp. (mm) 


3.0 


0.06 


2.8 


0.04 


2.6 


(T (mm) 


1.0 


0.0 


0.9 


0.0 


0.9 



Table 4. Regularity, Quality (sec): mesh repair times for both phases, in seconds; 
% nodes, nodes/2105: fraction and number of nodes moved by the repair pro- 
cedure; Mean, Max, a (mm): mean and maximal nodal displacements, standard 
deviation, in millimeters. 



Mesh quality distribution measured on the 1754 elements hip Atlas mesh 
is shown in Fig. |5]-a, and Fig. [B]-b gives the mean quality distribution in the 5 
generated patient-specific meshes, along with standard deviations. 

As in the previous illustrative case, the MMRep algorithm successfully worked 
for all 5 patients. All patient specific meshes were generated in less than 25 sec- 
onds and had submillimetric surface representation accuracy which remained 
unchanged by the repair procedures affecting less than 0.7% of the nodes. The 
maximal surface representation errors reported here are due to lack of refinement 
in the template mesh but also to the presence of noise in the digitized point 
clouds. Indeed, if the hand-held digitization pointer is lifted from the bony sur- 
face during the hip surface acquisition process, erroneous points can be recorded. 
These outliers are eventually averaged out during the elastic registration process, 
but their presence is revealed by the surface representation error measures. 

This second example demonstrates the adequacy of the MMRep technique 
in situations where only a fraction of the organ anatomy is known prior to 
modeling. The elastic registration process, thanks to an a priori knowledge about 
the organ of interest carried by the Atlas mesh, makes it possible to generate a 
FE model with high surface representation accuracy in the digitized zones and 
an approximate yet realistic organ shape in regions where no data is available. 
The repair phase ensures that the produced model meets the required quality 
standard and is suitable for FE analysis. 

4.3 Skin and bone face modeling 

This last use case demonstrates the application of the MMRep procedure to 
patient specific FE mesh generation in the context of orthognathic surgery, where 
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FE analysis helps predict the consequences of the intervention on the patient's 
features and facial expressions by simulating the effects of the repositioning of 
the jaw, maxillary or malar bones |37|38|8) . 

Mesh registration procedure In this application a manually assembled 3- 
layers mesh developed by Nazari et al. [39], shown in Fig. El is used to model the 
face muscles and fat. It is made of 8746 nodes forming 6030 hexahedrons and 314 
wedges. As we wish to model the interaction between face bones, muscles and 
features, the patient specific FE model must fit both the skin and skull recon- 
structed from each CT volume. To this end the Atlas outer and inner layer nodes 
are labeled "skin" and "bone" respectively and the multi-labels formulation of 
the elastic registration discussed in §2.71 is used. The inside Atlas mesh nodes, 
defining the inner layers of the face tissues, are unlabeled and follow the overall 
elastic deformation driven by the registration of the nodes labeled "skin" and 
"bone" . 

For the 50 patients included in this study (data provided by the MAP5 
Laboratory, University of Paris V), the bone and skin layers were segmented 
in the CT volumes using the Hounsfield scale and the resulting surfaces were 
reconstructed [40] and oriented along the anatomical axes in accordance with 
the generic face model. The Atlas mesh was aligned on the patient data using 
the nose tip position. Fig. [TU] shows sample patient skin and bone surfaces along 
with the translated Atlas model before elastic registration. 

The generic face mesh represents a subset of the complete patient head data, 
therefore, as discussed in the computed registration R is the one that fits 
the labeled Atlas mesh nodes to their corresponding destination skin or bone 
surfaces, producing the transformation that has to be applied to the generic 
mesh to specialize it for the specific patient. From the computational point of 
view, as a distinct set of distance maps has to be computed for each patient's 
skin and skull destination surfaces, the mesh generation process requires more 
time than in the previously described cases. Yet, in a pre-operative simulation 
scenario, computational delays are less critical than in an intraoperative FE 
analysis context. 

Results Unlike the previous use cases, this section does not give the detail 
of MMRep performance figures for each of the 50 patients. Instead, Table [5] 
summarizes the overall performance of our technique by presenting the mean 
registration speed, mesh repair cost and surface reconstruction accuracy. The 
surface representation errors shown here are the final measures performed after 
the repair procedure has been applied to the mesh. The nodal displacements 
amplitudes were evaluated separately for bone and skin layers. 

The impact of the mesh repair procedure is very limited as it only affects 
an average of 2% of the 1715 bone nodes and 0.3% of the 2180 skin nodes in 
the mesh. As for the surface representation, submillimetric mean accuracy is 
achieved for both layers. 
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Mean 


Max 




Elastic registration (sec) 


32 


96 


19 


Regularity (sec) 


28 


50 


12 


Quality (sec) 


4 


21 


2 


Bone 








Surface err. (mm) 


0.6 


13.8 


0.5 


Moved nodes / 1715 


39.8 


105 


23.9 


Displacements (mm) 


0.8 


3.8 


0.5 


Skin 








Surface err. (mm) 


0.4 


23.5 


0.4 


Moved nodes / 2180 


7.6 


45 


10.3 


Displacements (mm) 


0.3 


3.1 


0.2 



Table 5. Elastic mesh registration, regularization and quality optimization 
times, in seconds. For both bone and skin layers: final surface representation 
mean errors, number of nodes corrected by the repair procedure in each layer 
and nodal displacements amplitudes, in millimeters. 



The reported maximal errors are located in areas where the Atlas mesh lacks 
refinement. These regions clearly appear on Fig. 1111 where the surface represen- 
tation mean errors computed over the 50 cases are displayed as error maps on 
the initial Atlas bone and skin layers. The lightest areas represent mean errors 
between and 1 mm and the darkest areas errors above 3 mm. 

Fig. [TT] shows that the maximal skin layer errors are located around the ears 
which are clearly absent from the generic mesh, as can be seen in Fig. [9]-c. As 
for the bone layer, the maximal errors appear near the sphenoid bone and the 
zygomatic process as both regions have a very approximative representation in 
the Atlas mesh. 

Mesh quality distribution measured on the 6344 elements face Atlas mesh is 
shown in Fig. [I2\a. and Fig. [T2lb gives the mean quality distribution in the 50 
generated patient-specific meshes, along with standard deviations. 

In all cases, the MMRep algorithm was able to generate a patient specific 
bi-layer mesh suitable for FE analysis within a couple of minutes. The produced 
meshes exhibited submillimetric mean surface representation accuracy on both 
skin and bone layers, with larger errors localized around features absent or ill- 
defined in the generic mesh. 

Fig. [12] shows four examples of meshes fitted onto distinct patient mor- 
phologies. The Atlas face mesh used here was constructed based on a unique 
prognathicjfl patient, yet it was successfully used to model both prognathic and 
retrognathicu patients, which demonstrates the versatility of the registration 
procedure. 



® Having the jaws projecting beyond the upper part of the face. 
^ Having a mandible located posterior to its normal position. 
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Finally, Fig. [14] shows 25 thumbnails of registered face meshes demonstrating 
the variety of clinical cases embraced by the study. The two upper rows show 
retrognathic cases, the middle row average patients and the two lower rows 
prognathic morphologies. 

5 Discussion and conclusions 

A fast and automatic mesh generation procedure - the MMRep algorithm - has 
been presented, based upon the elastic registration of a generic, or "Atlas" , mesh 
towards patient specific structures, coupled with a robust mesh repair technique 
which ensures that element quality standards are met and FE analysis can safely 
be carried out on the resulting domain discretization. 

To our knowledge this is the first evaluation of a FE mesh registration tech- 
nique carried out on a wide range of 60 clinical data sets, illustrating 3 distinct 
use cases, raising both pre- and intraoperative biomechanical modeling issues 
and relying on fully or partially available patient data. In all situations the MM- 
Rep technique automatically generated a patient specific quality mesh within 
minutes, which strongly contrasts with time-consuming manual mesh assembly 
procedures involving a human expert operator. 

In all studied cases the regularity of the elastic deformation preserved the 
mesh elements from excessive distortions and the repair algorithm was able to 
find a suitable nodal configuration while maintaining a satisfactory surface rep- 
resentation accuracy. Only a small fraction, less than 1%, of the mesh nodes 
positions needed to be corrected by submillimetric displacements. 

Furthermore, the elastic registration formulation made it possible to compute 
a mesh deformation driven by multiple anatomical structures or sub-structures, 
as shown in ^4^^ This feature could easily be used in the femur model generation 
scenarii, §4.1l and l4.21 for example for distinction between cortical and spongious 
bone layers, which have distinct mechanical properties. 

This study shows that the proposed elastic registration technique is well- 
suited to the addressed meshing problem: computational times are short and 
submillimetric surface representation accuracy is achieved. However the MMRep 
procedure could be supported by any elastic registration algorithm |41|42|27] . 
provided that the following properties are ensured: 

— A smooth deformation field, ideally C^, can be estimated by the registration 
algorithm within reasonable computational times. 

— Non-folding and bijection of the registration are ensured. 

— Accurate registration inverse can be computed if the patient data is regis- 
tered onto the Atlas and the inverse deformation is applied to the generic 
mesh in order to make it specific, as described in ii4.1l and 14.21 

— The simultaneous registration of different structures, if required by the biome- 
chanical modeling, can be computed as shown in §4.31 

The Atlas based approach presented here relies on the definition of a generic 
model of the target organ with the desired elements layout and possibly some 
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identified sub-structures of interest. Tliis modeiing effort oniy needs to be done 
once but the resuiting Atias mesli must be carefuiiy designed so tiiat tlie ana- 
tomicai features of interest can be properly registered to tlieir patient specific 
counterparts. 

The tests carried out on our database suggest algorithm robustness although 
no formal proof of convergence could be given. In extreme cases where strongly 
distorted pathological organs diverge from the average Atlas shape, mesh regu- 
larity can be lost and the repair strategies proposed here may fail. This issue can 
be overcome by working with a pool of Atlas meshes that represent the main 
deformation classes likely to be encountered. Such Atlas variations can be easily 
generated by successfully applying the MMRep procedure to a representative 
case and using the resulting quality mesh as a starting point for subsequent 
mesh registrations of similar configurations. 

A number of enhancements to the presented MMRep technique can be fore- 
seen. Breaking the sequential registration-repair scheme, the algorithm could 
benefit from the integration of the repair process within the elastic registra- 
tion computation itself. This could be implemented either as a new energy term 
replacing the ad-hoc mechanical formulation controlling the deformation regula- 
rity, or as a per-iteration post-processing callback which, although more straight- 
forward, would unfortunately raise the issue of the convergence of the elastic 
registration algorithm. 

Future works also include the broadening of our mesh repair approach by 
taking into account other element types such as pyramids and tetrahedrons, as 
well as other quality criteria, such as the face warping factor I43j that measures 
each element's face nodes coplanarity. Stronger constraints on the evenness of 
the generated meshes should enhance the numerical stability of the FE analysis 
carried out. 

Finally, we can imagine an ideal FE mesh generation algorithm that performs 
mesh registration and repair directly in patient 3D image space. This fully in- 
tegrated organ modeling tool could be achieved by replacing the distance based 
elastic registration formulation used here by an appropriate segmentation energy 
suitable for the considered imaging modality. 
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(c) (d) 

Fig. 5. Sample patient specific mesh. Proximal epiphysis: (a) segmented patient 
data (black dots) and Atlas mesh; (b) patient specific mesh fitting the femur 
head surface. Distal epiphysis: (c) segmented patient data (black dots) and Atlas 
mesh; (d) patient specific mesh fitting the condyles. 
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Fig. 6. Mesh quality statistics for Atlas (a) and patient-specific (b) full femur 
models. 




Fig. 7. Intraopcrativcly acquired point clouds for the 5 considered patients 
(black dots). To localize the intraoperatively accessible bone region, a surface 
mesh of the femoral head is showed with each data set (wireframe surface). 
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Fig. 8. Mesh quality statistics for Atlas (a) and patient-specific (b) hip models. 
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Fig. 10. (a) Skin and skull surfaces reconstruction from a sample patient CT 
volume, (b) Generic face mesh rigidly registered with the patient data. 
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Fig. 11. Face mesh registration mean errors represented as color maps, (a) Error 
maps color code in millimeters. Skin layer: (b) anterior view; (c) lateral view. 
Bone layer: (d) anterior view; (e) lateral view. 
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Fig. 12. Mesh quality statistics for Atlas (a) and patient-specific (b) face models. 
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Fig. 13. Four examples of the generated face raodels. For clarity only the skin 
surface segmented from the CT volume is shown here although the produced FE 
models also fit to the underlying skull surface. 
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Fig. 14. Sample of 25 registered face meshes. For each patient, the transparent 
skin surface mesh is shown along with the registered FE mesh. For clarity the 
skulls were omitted. 



